Onset of fiuidization in vertically shaken granular material 



o 

Q 



Thorsten Poschel 1,2 , Thomas Schwager 2 , and Clara Saluena 2 
1 ICA Universitat Stuttgart, Pfaffenwaldring 27, D-70569 Stuttgart, Germany 
2 Humboldt- Universitat zu Berlin, Institut fur Physik, Invahdenstrafie 110, D-10115 Berlin, Germany 

(February 1, 2008) 

When granular material is shaken vertically one observes convection, surface fiuidization, spon- 
taneous heap formation and other effects. There is a controversial discussion in literature whether 
there exists a threshold for the Froude number F — AoU>o/g below which these effects cannot be 
observed anymore. By means of theoretical analysis and computer simulation we find that there 
is no such single threshold. Instead we propose a modified criterion which coincides with critical 
Froude number r c = 1 for small driving frequency u>o- 

PACS: 81.05.Rm, 83.70.Fn, 46.30.My 
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I. INTRODUCTION 

When granular material in a rectangular container 
is imposed to vertical oscillation under certain condi- 
tions one observes a variety of effects, such as convec- 




spontaneous heap 
llLOl, oscillons flill 



tion jj]-[3), surface fiuidization 
formation |} 10 1, surface patterns 
and others. The common feature of all these effects is 
that particles change their position with respect to each 
other. Provided the particles do not change their me- 
chanical properties during observation time (by polish- 
ing, comminution etc.) the condition for this motion is 
that neighboring particles separate from each other at 
least for a small part of the oscillation cycle T = 2-k j w®. 

There is a controversial discussion in the literature 
whether there is a critical value of Froude number 



r, 



(1) 



below which the above mentioned effects vanish, with 
Aq and loq being the parameters of the sinusoidal mo- 
tion of the container. In many experimental observa- 
tions (e.g. |0,^,0Jl^J^^,^|) and computer simulations 
(e.g. |15| , |l6| ) such a critical number T c was found. Several 
authors believe that the value is T c = 1. In numerical 
simulations, however, surface fiuidization and convection 
has been found for T < 1 |p|,p|Jl7|. Therefore, some au- 
thors believe that T is not the proper criterion to deter- 
mine the degree of fiuidization of a granular system . 

In this article we discuss the response of granular ma- 
terial with respect to vertical oscillation in the limit of 
an one dimensional approach: the lowest bead of a verti- 
cal column of N identical spherical beads is shaken with 
zq = Aq cos uot and the other beads move due to their in- 
teraction force and gravity g. We study the motion of the 
entire column and can show that particles can lose con- 
tact to their neighbors even for the case that V = A uiq / g 
is significantly less than 1. 

Adjacent spheres k and k + 1 of radius r and mass m 
at vertical positions Zk and Zk+\ interact with their next 
neighbors by 



F, 



k,k+l = 



t 3/2 



fc+i 



(2) 



with fj, and a being elastic and dissipative material con- 
stants, i.e. functions of Young modulus, Poisson ratio 
and dissipation rate (for details of the derivation of Eq. 
(||) see [|9|). £ is the compression 2r — \zk — %k+i\ of 
the spheres. The height of the column is L = 2Nr. Ex- 
pression (J2J) is valid if the typical relative velocities of 
adjacent spheres are far below the speed of sound in the 
material of the spheres. Certainly this condition holds 
for typical vibration experiments. 

Introducing new coordinates Uk = Zk — 2rk (k = 
. . . N) the compression of two adjacent spheres is 



£,k,k+l — Uk — Uk+1 



(3) 



Applying these definitions in Eq. (g) and adding gravity 
g we get 



1 



Zk 



Fi 



k,k+l 



~n 

-/j,^/r(u k - ufc+i) 3/2 



(4) 



fr (ii k - Uk+i) y/u k - u k +i 



The 0-th sphere is fixed at the oscillating table, hence its 
position is 

Zo{t) = Uo(t) = AoCOSLUot. 

We are interested in the critical parameters of driving 
(Aq, wo) when the JV-th particle loses contact, i.e. when 
un > Un-i- We define the "response" i?(wo) as the ratio 
An/Aq where An is the amplitude of the iV-th particle 
at frequency loq and Aq is the amplitude of the driving 
vibration. R(ljq) can be calculated by convoluting the 
motion 2jv(t) with exp(iuJot) . Suppose A^oJ^/g > 1 the 
iV-th particle separates from the N — 1-st. If we would 
find Aq < Apf the critical Froude number T c — AQLOg/g 
would be less than 1. We will show that there is a range 
for uq where this is the case. 
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In the next section we will formulate the problem in a 
continuum approach and derive a nonlinear partial differ- 
ential equation for the motion of the column of particles. 



This equation is solved in section [II in the limit of elas- 



tic material properties, i.e. by dropping the dissipative 
terms. Once the solution for the elastic case has been 
discussed in detail, it is easier to study the influence of 
the dissipative term and to derive the solution of the full 
equation of motion which is done in section |y|. Section 
^ compares the analytic results with a molecular dynam- 
ics simulation of the original (discrete) problem stated in 
Eq. (0). Finally, we discuss the results. 



II. CONTINUUM APPROACH 

To study the system analytically we use an one dimen- 
sional continuum approach. To this end we perform a 
Taylor expansion of the force with respect to the radius 
r and consequently consider the limit r — > 0, N — > oo 
with 2rN = L = const. First we have to replace the dis- 
placements Uk by u(2kr) introducing the displacement 
field u(z) which is a continuous function of z. With Eq. 
(||) we find from Taylor expansion 

£ fe ,fc+i = u{2kr) - u(2kr + 2r) = - 2r [du/ dz)\ z=2kr . 

The net force experienced by the fc-th particle is 

Fk = Fk^k+l — Fk-l.k 



—ay/r \ ^k,k+i\/£,k,k+i - €k-i.ky/£k-i,k 



-2V2r 2 n 



du k 
dz 



.3/2 



+2v / 2r 2 a 



d 2 u k 
dtdz 

with the abbreviations 



duk 
dz 



dz 

d 2 u k ~ 
dtdz 



3/2' 



du k - 



Uk = u(2kr) 
duk du 
dz dz 



(5) 
(6) 



2k r 



Both expressions in square brackets are expanded again 
and Eq. (§) turns into 



m 



3\/2 



du 
dz 



3/2 



With 



3\/2u 



dtdz 



3V2c 
irp 



du 
dz 



the continuum formulation of 



dt 2 9 dz 



du 

dz 



_du\ 
~~dz~) 



3/2 



dtdz 



du 
dz 



(7) 



z=L 



where g accounts for the gravitational force. 

III. LIMIT OF ELASTIC PARTICLES 

In the following we consider Eq. ([?]) the limit of no 
damping ((3 = 0). Using new variables 



x = 1 . 

L 

7-5 \ 1/6 



1/6 



gn 



Eq. (Q) turns into 



d 2 u 2 ^1 3 

<9t 2 ^ 7 dx 



du\ 3/2 
dx) 



du 
dx 



= o. 



(8) 
(9) 



Equation (B) is defined in the range x S [0, 1]. The time 
independent solution U(x) of (H) is 



U{x) = ^ 7 2 (^ 5/3 - l) 



(10) 



The solution of (|§) can be considered as a superposition 
of the static solution ( |l0[) and a perturbation w(x, r). 
Inserting u = U + w in dq) we find: 



d 2 w 2 1 d 

dr 2 ^ 7 dx 

2 1 d 

7 dx 



dll dw 
dx dx 

dU \ 3/2 

dx 



3/2 

3 



/ dU dw 
dx dx 



3d_ 

2dx 



i dw 



dx 



(11) 



By separation of variables w = T(t,SY)X(x,{Y), i.e. a 
standing wave Ansatz we obtain two ordinary differen- 
tial equations for T and x: 



T = -n 2 T 



3 d 



-1/3 



dX 



= -n 2 x. 



2 dx \ dx 
with fi being a real number. For T(r, f2) one gets 
T ~ exp(ifir) . 



(12) 
(13) 
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The solution of the spatial equation ([13]) can be found 
using the Ansatz 

X(x,n) = xV 3 f(y), y=^V6nx 5 / 6 , 



which yields 



dy 2 dy 



(14) 



Eq. (14) is the Bessel equation of order 2/5. Hence the 
solution of (O) is 



X(x,Q) 



1/5 



r 1 1 1 n 2 /V/ 3 j_ 



(15) 



An expression containing J 2 / 5 would be a solution too, 
however, it does not satisfy the condition (^|). The pref- 
actor in ( |l5| ) has been chosen to assure A(0, O) = 1. 
Hence the solution for a single vibrational mode mq is 



(16) 



mo = exp (iftr) A(x, f2) . 



Without prior knowledge the full solution of Eq. (Ill 
has to be assumed to be a superposition of vibrational 
modes for all real (rescaled) frequencies ft: 



dft,4(ft)exp (iftr)X(:E,ft) . 



(17) 



In the steady state of pure sinusoidal excitation of the 
base, i.e., when all non-oscillatory perturbations which 
originate from the initialization have been damped out, 
Eq. (O) is the full (steady state-) solution of Eq. (|TTj) . 

The function A(ft) represents the excitation of the 
mode at frequency ft. The boundary condition at the 
top of the chain is automatically satisfied, whereas the 
boundary condition at the bottom reads 



(18) 
(19) 



u(1,t)= dO J 4(n)exp(ifir)X(l,0) 

J — oo 
= A COS FIqT . 



One can see that the integrand of Eq. ( |l8| ) can be nonzero 
only for ft ^ fto. This means that for ft ^ fto either A(ft) 
or X(l,Sl) have to be zero, i.e. for all frequencies for 
which X(l,ft) is nonzero the respective amplitude must 
be zero, whereas for all frequencies which are a root of 
X(1,Q,) = the amplitude can be nonzero. Therefore, 
we find that the full solution of Eq. (jll]) is a superposi- 
tion of the vibrational mode of the frequency of shaking 
Oo and of a discrete set of frequencies ft*, (k = 1 . . . oo). 

Note that ftfc are no rational multiples of each others 
since the roots of Bessel-functions are incommensurable 
(see Eq. fll5])). Therefore, to determine the maximum 



acceleration of the topmost particle it is sufficient to con- 
sider only the mode of the external excitation. All other 
vibrational modes can only further increase the maximal 
acceleration. 

The above defined response R is the ratio A^/Aq. 
Since the zeroth particle corresponds to x = 1 and the 
N-th to x = we can write 



ir x (ft ) - 



x(i,n ) 



X(0,Qo) 



= \x(i,n )\ 



1/5 



ft 



2/5 



(20) 



The response R is an amplification factor, hence, the 
value g/R(£lo) is the critical acceleration of the driving 
vibration pC)| . R is larger than 1 for all driving frequen- 
cies <x>o . This means that for any driving frequency u>q and 
driving amplitude A the amplitude of the top particle 
of the column Ajq at frequency ujq will be larger than A . 
Therefore, for A^uy^jg = 1, i.e. when the 7V-th particle 
separates from the N — 1-st, we find Aouj^/g = T c < 1. 

According to the above arguments we have to replace 
the condition T > 1 which was supposed to be the pre- 
condition for surface fluidization, convection etc., by 



A ul/g = T>R- 1 (lo ) ■ 



(21) 



The function R 1 (wo) over ujq is drawn in Fig. g (dash- 
dotted line, -RJ/ 1 )- For the system parameters we used 
Aq = 0.01 mm, elastic constant k ~ 2. 8 TO 4 m 2 /sec 2 (rub- 
ber with Young modulus Y = 4- 10 7 Pa) and L — 0.6 m. 
The curve reveals pronounced resonances at Eigenfre- 
quencies LOk where R~ l becomes minimal (only the first 
resonance is shown in Fig. Q). 

All experiments on surface fluidization and convection 
which can be found in literature have been performed far 
below the first resonance. Therefore, of particular inter- 
est to practical purposes is the limit of small frequency 
wo, i.e. below the first Eigenfrequency. The Taylor ex- 
pansion of Eq. (p0|) yields R^ 1 (flo) for small Qq 



R- 1 = 1 - -ft 2 



I-* 

5 



1/3 



(22) 



Given the container vibrates with frequency luq. Then 
for the critical amplitude Aq of the vibration when the 
top particle separates, i.e. when the material starts to 
fluidize, one finds 



A = 



gx? 



1/3 



(23) 
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FIG. 1. The response function R^ 1 due to Eq. ( p(i| ) for 
a column of elastic particles (dash-dotted). For dissipative 
particles (/3 =127m 2 /sec): circles, R^ : numerical integra- 
tion of Eq. M) at small amplitude; full line, R^ co : analytical 
solution (po|) of the full Eq. (Q) including dissipation; dashed 
line, F c (o;o): result of a direct simulation of Eq. ([|) (expla- 
nation see text) 

Surprisingly even for very small frequencies where 
R~ x —> 1 one finds that the critical amplitude is reduced 
by a constant as compared with g/uj 2 ,. So although the 
value of the response function comes arbitrarily close to 
one, the critical amplitude differs from the expected one 
by a constant. However, this does not mean that the 
critical Froude number becomes a constant. 

Eq. (|l^) describes the behavior of the column of grains 
for the case of purely elastic contact (a — 0). If the dis- 
sipative material properties are taken into consideration 
the full equation (0) has to be solved which will be dis- 
cussed in the following section. 



IV. DISSIPATIVE PARTICLE INTERACTION 

We will consider, as before, small perturbations w 
about the static deformation of the chain under gravity 
which propagate from the bottom. The dissipative term 
is characterized by the parameter /3 in Eq. (0). From this 
equation, again introducing the static solution given by 
Eq. ( |Io|) and by using the same transformation for the 
spatial coordinate x = 1 — z/L, one obtains the corre- 
sponding linearized wave equation for dissipative mate- 
rials, 



d 2 w 



( 9 y 

\kL 5 ) 



3 d_ 

2 K dx 



13 ■ 



. dw 

dx 



d_ 

dx 



X 3 



d 2 w \ 
dtdx) 



(24) 



In this case it is less useful to introduce the rescaled time 
variable r, while it proves convenient to define 



nil- 
9 



5\ 1/6 



(25) 



being the natural length scale coming out from the anal- 
ysis. By means of the Fourier transform, 



W(x,u) 



1 



e- luJt w{x,t)dt, (26) 



equation (g4() becomes 

f/t-iw/3 d ( 1/3 dW 



■WW 



which has the same structure as Eq. (|l3|). Hence the 
same transformations apply in this case and the general 
solution reads finally 



W(x,u>) = x 



hi T 5/ 6 



|k — icu/3 



C 2 J 



(28) 



The part of the solution depending on J 2 /5 carries a di- 
vergence at x = (z = L), and therefore Ci = is 
required for the solution to be physical. The condition of 
the free end at x — (z — L) is satisfied automatically, 
as in the case j3 = 0. The solution has exactly the same 
structure as the solution of the elastic problem, Eq. (jig) , 
and the only change is that the argument of the Bessel 
function has an imaginary part. If one considers only 
the mode luo, which propagates from the bottom [x = 1) 
with amplitude A$, the solution reads 



w(x,t) = Re < 



J_ 



A e iuot x* 



|k — iuo/3 



r 5/6 



(29) 



The fluidization condition at the top of the chain 
\d 2 w/dt 2 (x — 0,t)\ > g can be written in general 



< 



5 

u 2 
9 



uj e 



2/5 



uj £ 



(30) 



4 



Since the Bessel functions of the first class only have 
zeros on the real axis, i? _1 (a;o) can no longer be zero 
for any frequency if (3 ^ 0. This means that the sharp 
resonances displayed by R(u>q) when (3 = disappear 
and are replaced by more or less pronounced minima 
in dependence on the damping constant (3. This can 
also be observed in Figure [l] (full line). Increasing val- 
ues of (3 make the response smoother, deviating from 
that of the elastic case earlier, and the local minima 
translate on the frequency axis appreciably, see Figure |^. 



DC 




50 

co [l/sec] 

FIG. 2. The response function over the frequency for dif- 
ferent dissipative constants: f3 = (0, 100, 200, . . . , 1000) 
m 2 /sec (bottom to top). With increasing damping the min- 
imum becomes less pronounced and shifts to lower frequen- 
cies. 



Analogue to the elastic case (Eq. (|22|)), for small fre- 
quencies R~ 1 (u>q) can be expanded into a Taylor series: 



5 \gn 
L 5 \ 2/3 

gn 2 



1/3 

3 

100 



4--) 



1/3 



(31) 
(32) 



The contribution due to the dissipative parameter en- 
ters the Taylor expansion at the fourth power of the fre- 
quency. Therefore, the analysis of the elastic case given 
by Eq. (^||) remains valid for small frequencies. 

It is interesting to note that due to Eq. ( |3l| ) 
there always exists a global minimum below the value 
i? _1 (o;o) = 1, regardless of the value of the dissipa- 
tive parameter. We want to study for which range of 
frequencies the inverse response is smaller than 

one for varying damping constant (3: The lower bound- 
ary of this interval obviously is w = 0. To deter- 
mine the upper boundary w max we solved numerically 
the equation i? _1 (aj max ) = 1 for different values of 
(3. The result of this calculation is shown in Fig. [|. 




200 400 600 800 

P [m /sec] 

FIG. 3. The frequency oj max at which the inverse response 
function becomes larger than one is almost constant with 
varying dissipative parameter f3 as long as the damping is 
not to small. For very low dissipative parameter j3 one finds 
a non-smooth function (small figure). 

One can see that for high enough damping this fre- 
quency w max varies only slowly with (3. The curve al- 
most saturates at w max w 200 sec -1 which is close to 
the first maximum of the undamped inverse response 
(see Fig. ^, dash-dotted line). Although the valley 
of the inverse response function becomes smaller with 
increasing damping (Fig. |3|), even for larger damping 
_R _1 is smaller than one in a finite frequency interval, 
i.e. the effect of amplitude amplification exists for al- 
most the entire range of frequencies between zero and 
the first maximum of the undamped inverse response. 




400 

oo [l/sec] 

FIG. 4. The response function for different values of the 
dissipative constant (3 =10, 20, 30,. . . ,100 m 2 /sec (bottom 
to top) together with the elastic curve (3 = (dash-dotted). 
With decreasing j3 the curves are influenced by higher order 
minima of the response function. This explains the steps in 
the curve drawn in Fig. ^| for small values of the dissipation 
(3. 
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For small enough damping the frequency range of am- 
plitude amplification R" 1 (uj ) < 1 extends beyond the 
position of the first maximum of the undamped inverse 
response. Figure ^ shows the response function for dif- 
ferent dissipative parameter (3 together with the elastic 
case ((3 = 0, dash-dotted). The frequency wo at which 
i? _1 (/3) = 1 increases with decreasing damping. In this 
Fig. ^ one can clearly see the extension of the range 
of amplitude amplification. If we decrease the damping 
parameter, starting at high values where the amplifying 
range is limited to the first "well" we cannot expect a 
large change since, as long as we remain limited to the 
first well, there is a upper bound (the frequency of the 
first maximum of i? _1 ) to this range. Reducing j3 further 
we will eventually reach values for which the amplifying 
range extends to the second well. Even for values of 
which are only slightly below this threshold the range 
will now span again almost the entire range of the sec- 
ond well up to the frequency of the second maximum. So 
there will be rather sharp steps in the dependence of the 
upper range limit of the damping instead of only gradual 
changes. This behaviour explains the steps in the closeup 
of Fig. |. 



the column N can separate from the N — 1-st even if the 
container is oscillated with Aquiq/q < 1. As main result 
we derived a modified condition for the topmost particle 
to separate from its neighbor. We could show that in- 
stead of the widely accepted condition T c = Aqloq / g > 1 
one has to satisfy AquJq/ 1 g > R^ 1 where is a function 
of ujq. We have shown that independent on the mate- 
rial properties there exists always a range oj E [0. 
for which the amplitude of vibration Aq is amplified, i.e. 
for which the top particle can separate (the material flu- 
idizes) even if Aqluq j g < 1 . Numerical calculations agree 
well with the analytic results. 

Whereas the critical Froude number r c > 1 is cer- 
tainly the proper criterion to predict whether a single 
rigid particle will jump on a vibrating table, we suspect 
that this number is not suited to describe the behavior 
of a vibrated column of spheres, and the more it is not 
a criterion for surface fluidization of a three dimensional 
granular material. 

The authors wish to thank E. Clement, N. Gray, 
H. J. Herrmann, H. M. Jaeger, S. Luding, S. Roux and 
L. Schimansky-Geier for helpful discussion. 



V. NUMERICAL RESULTS 



To check the analytical results and in particular the va- 
lidity of the continuum approach we calculated for a 
finite value of damping (3 from the numerical simulation 
of Eq. ([|). The circles in Fig. [I] display the reciprocal re- 
sponse over uq with fixed amplitude Aq = 0.01 mm, 
elastic constant n — 2.8 • 10 4 m 2 /sec 2 (rubber with Young 
modulus Y = 4 • 10 7 Pa) and L = 0.6 m. Fig. [l] shows 
that for small frequency loq and small damping a the un- 
damped theoretical curve (dash-dotted) agrees well with 
numerical data. If one compares the numerical result 
with the damped solution according to Eq. (|30| ) (full line 
in Fig. [j]) the agreement with theory is very well. 

To check the validity of our linear theory we also de- 
termined directly by integrating Eq. (^) at which Froude 
number T c the particles start to jump. The results of this 
calculation are shown in the dashed curve in Fig. and 
agree well with the linear theory. To obtain the value 
of r c for a given frequency loq we determined an upper 
bound A^ for the critical amplitude where one observes 
jumping and a lower bound Aq where no jumping occurs. 
Then we narrowed the interval Aq — Aq by testing an am- 
plitude between Afi and Aq until (A^-Aq )/A£ < 10~ 3 . 

VI. DISCUSSION 

For the case of a vertical column of viscoelastic spheres 
we derived a linear wave equation in one dimensional ap- 
proximation. We have shown that the sphere on top of 
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